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Abstract. We use the usual method of Schwarzschild to construct self-consistent solutions for the triaxial de Zeeuw & Carollo (1996) 
models with central density cusps. ZC96 models are triaxial generalisations of spherical y-models of Dehnen whose densities vary as 
r"'' near the center and r^* at large radii and hence, possess a central density core for y = and cusps for y > 0. We consider four 
triaxial models from ZC96, two prolate triaxials: (p, q) = (0.65, 0.60) with y = 1.0 and 1.5, and two oblate triaxials: (/?, q) = (0.95, 0.60) 
with 7 = 1.0 and 1.5. We compute 4500 orbits in each model for time periods of lO^To- We find that a large fraction of the orbits in 
each model are stochastic by means of their nonzero Liapunov exponents. The stochastic orbits in each model can sustain regular 
shapes for ~ lO^T/j or longer, which suggests that they difi"use slowly through their allowed phase-space. Except for the oblate triaxial 
models with y = 1.0, our attempts to construct self-consistent solutions employing only the regular orbits fail for the remaining three 
models. However, the self-consistent solutions are found to exist for all models when the stochastic and regular orbits are treated in the 
same way because the mixing-time, ~ 10*Td, is shorter than the integration time, lO^T/j. Moreover, the "fully-mixed" solutions can 
also be constructed for all models when the stochastic orbits are fully mixed at 15 lowest energy shells. Thus, we conclude that the 
self-consistent solutions exist for our selected prolate and oblate triaxial models with y = 1.0 and 1.5. 
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It is natural to assume that elliptical galaxies can be triaxial, 
such as ellipsoids (Binney 1978). In this regard, Schwarzschild 
(1979) developed a method to explore the existence of triax- 
ial ellipticals using a catalogue of numerically integrated orbits, 
which eventually became a popular tool. The triaxial models con- 
sidered by him feature a central density core, and also show 
isophotal twists (Chakraborty & Thakur 2000) representing the ob- 
served photometric properties of the majority of elliptical galaxies 
(Jedrzejewski 1987; Peletier et al. 1990). However, the ground- 
based (Moller, StiavelU & Zeilinger 1995) and Hubble Space 
Telescope (Crane et al. 1993; Jaffe et al. 1994; Ferrarese et al. 
1994; Lauer et al. 1995; Faber et al. 1997) observations reveal 
that instead of central density cores, most of the elliptical galaxies 
have central density cusps. Thus, Merritt & Fridman (1996) stud- 
ied the effect of central density cusps on triaxial configurations but 
their models do not show isophotal twists. In their study, it was 
found that a "fully-mixed" solution exists for a maximally triaxial 
model with a weak-cusp. The aim of this paper is to construct self- 
consistent solutions for more realistic triaxial models in the sense 
of representing both of the above-mentioned significant observed 



properties of elliptical galaxies, namely central density cusps and 
isophotal twists. Here we used the triaxial models given by de 
Zeeuw & Carollo (1996) (hereafter ZC96), which have a central 
density core for y = and a cusp for y > 0, and also show isopho- 
tal twists in their projection along the line-of-sight. Furthermore, 
it is worth to mention that the ZC96 triaxial models are found to 
be useful in constraining the intrinsic shapes of elliptical galaxies 
using their projected properties (Thakur & Chakraborty 2001). 

The remainder of this paper is organized as follows. In Sect. 
2, we present the ZC96 triaxial models. Sect. 3 describes equa- 
tions of motion for the ZC96 triaxial models. Sect. 4 presents the 
integration of orbits. Sect. 5 deals with the computation of the 
Liapunov exponents. The method to construct the self-consistent 
triaxial models, is presented in Sect. 6. Sect. 7 is devoted for the 
results and discussion. 

2. Triaxial de Zeeuw-Carollo models 

We have used the triaxial potentials of ZC96, given as 
V(r,6,(f>) = u(r) - v(r)Y°(e) + w(r)Y^(e,(f>), 



V(x,y,z) = M(r) - v(r)(2z^ 



■/)/2r2-fw(r)3(x2-/)/^2,(l) 
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where (r, 0, (p) are spherical coordinates defined such that x - 
r sin cos ^,y - r sin sin (p and z-r cos 9, the functions (^) - 
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I cos^ ^ ~ 2 and 0) = 3 sin^ 6 cos 20 are usual spherical har- 
monics, and u(r), v(r) and w{r) are three radial functions. Here u{r) 
is chosen to be the potential of the spherical y-models of Dehnen 
(1993), defined by 

™lnj^, for 7 = 2, 

(fe-^^)'-^-!], for 7^2, 

where M is the total mass of the model, is the scale-length 
and cusp parameter y can have a value in between < y < 3. 
Furthermore, the functions v{r) and w(r) are considered as follows: 



u{r) 



v(r) = - 



w(r) = - 



GMr^r^-y 
(r + r4)4-r' 



(3) 



(r + r2)4-r' 

where ri, ra, are constants. 

The associated density distribution p{r,0,(/>) follows from 
Poisson's equation 

pir, 6, 4>) = fir) - gir)Y"2(0) + h(r)Y^(e, 0), 

Pix, y, z) = fir) - gir)i2z' - - y^)/2r' + hir)3ix' - y^)/r\ (4) 

where fir), gir), and /i(r) are taken from Eq. (2.5) of ZC96. The 
four ratios ri /r,,, r4/ro can be expressed in terms of y, and of the 
axial ratios of the density distribution at small and at large radii, 
respectively denoted by ipo,qo) and (/7oo,<?oo), where the surfaces 
of constant density are approximately elhpsoidal, i.e., p ~ pim^) 
with = +y^ Ip^ + z^lq^. ZC96 models have a central density 
core for y = 0, while they have cusps in which the density diverges 
as r at small radii for 7 > 0. On the other hand, the density fall 
off as r"* at large radii. 

In this paper, we fix the values of axial ratios as Po = Poo = p 
and qo - Qoo = q for all the calculations. 

3. Equation of motion for triaxial de Zeeuw-Carollo 
models 



Equation of motion is given by 



dfl 



= F = -^Vix,y,z) 



(5) 



where F is the force per unit mass, and Vix,y,z) represents the 
triaxial potentials of ZC96 given in Eq. (1). 

The force components (F ^, F^, F^) = iFi,F2, F^) representing 
the equations of motion by three scalar equations in the Cartesian 
coordinates, needed at each step of an orbit integration, can be 
calculated by the partial derivatives of Vix,y,z) with respect to 
ix, y, z) = ix\ , X2, X3). These can be written as 



Xi 

Fi = -- 

r 



Iwd 2Wr\. 2 2^ A "^r 



(6) 

where (Ai , A2, A3) are equal to (1 , 3, 6), (1, 3, -6), and (-2, 3, 0) for 
i = 1, 2, 3, respectively. Furthermore, 



Ud 



l-y 



1 



vd = -n 



Wd 



r+l (r -I- 1)2 

(2 - y)r^-T (4 - y)^--^ 



ir + r2)^-y ir + r2f-y 
i2-y)r^-y iA-y)r^-y 



ir + rA)^-y ir + rA^-y 



(7) 



are the derivatives of Ur, Vr, and Wr with respect to r, respectively. 
Here the expressions of Uy, Vy, and Wr can be derived from those of 
M(r), v(r), and w(r) given in Eqs. (2)-(3) using G = M = rg = I, 
respectively. 

4. Integration of orbits 

To integrate orbits, we have considered four triaxial models of 
ZC96, which are listed in Table 1 . Here the first two models are 
prolate triaxials: ip,q) - (0.65,0.60) with 7 = 1.0 and 1.5, while 
the remaining two are oblate triaxials: ip,q) - (0.95,0.60) with 
7 = 1.0 and 1.5. In this paper, the former models will be called 
Models FT, whereas the latter models will be referred to as Models 
OT. 

Since the triaxial mass models of ZC96 are very centrally 
concentrated, the numerical algorithm required for integrating the 
orbits must be extremely accurate and flexible. We have used a 
FORTRAN routine, RK78, which was kindly made available by 
Prof. D. Pfenniger. This routine follows tiie 7/8 order Runga-Kutta 
algorithm described in Fehlberg (1968), which incorporates a vari- 
able time step in order to maintain a specified accuracy from one 
integration step to the next. The accuracy parameter REPS in all 
the integrations is chosen to be 10"**. Energy is typically conserved 
to a few parts in 10^ over 10^ dynamical times with this choice of 
REPS. 

For each of the four selected triaxial models, we have followed 
the scheme developed by Schwarzschild (1993) in assigning initial 
conditions from the x - z start-space with = = and sta- 
tionary (equipotential) start-space with Vx - v, = = 0. As in 
Merritt & Fridman (1996), the orbits in both start-spaces are as- 
signed a value from a set of 20 energies, defined as the values of 
the potential on the x-axis of a set of 20 ellipsoidal shells. The ra- 
dius, energy and dynamical time (Jc) of each ellipsoidal shell are 
given in Table 2 for all models. Here energy dependent "dynamical 
time" To is defined as the period of the 1 : 1 resonant x-y periodic 
orbit. For the x - z start-space, a total of 150 starting points per 
shell are calculated for all models. Moreover, the stationary start- 
space grid is defined as in Schwarzschild (1993). But, as opposed 
to 64 starting points in Merritt & Fridman (1996), only 25 starting 
points are calculated in each of the three sectors on an equipotential 
octant, resulting a total of 75 starting points per shell in the station- 
ary start-space. Here, as the Liapunov exponents are computed for 
long periods of lO^r^ (see Sect. 5), we reduce the number of start- 
ing points in order to save computing time. Thus, a total of 4500 
starting points are calculated for each of the selected models due 
to their division into a set of 20 elhpsoidal sheUs. For each of these 
starting points, we have integrated the orbit over a time interval of 
lO^Jz). 

5. Computation of Liapunov exponents 

For detecting and quantifying stochasticity, we have followed 
Merritt & Fridman (1996) and computed approximations to the 
six Liapunov exponents, which can be ordered by size as cr\ > 
(T2 > ... > CTg, by integrating each of 4500 orbits for long pe- 
riods of lO^To using the Gram-Schmidt orthogonalization tech- 
nique described by Benettin et al. (1980). In order to carry out 
this technique, we have used a FORTRAN routine, LIAMAG, 
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which was kindly made available by Prof. D. Pfenniger of Geneva 
Observatory group. In this routine, the second derivatives of the 
potential with respect to position are required for determination 
of the evolution of the perturbed orbits (Udry & Pfenniger 1988), 
which are given in Appendix A. 

As in Merritt & Fridman (1996), we have restricted our atten- 
tion to the three positive Liapunov exponents, since cr, = -o"7_, 
with i = 1, 2, 3. For all models, we have calculated the Liapunov 
times (tl), corresponding to the instabiUty timescale between 
neighboring orbits, which are found to be approximately same for 
all energy shells when they are scaled in dynamical times (To)- 
Table 1 gives these calculated values of Liapunov time (tl) in 
units of To for all models. Here Models PT with y = 1.5 ap- 
pears as the most stochastic one, since its Liapunov time (tl) is 
smaller compared to other selected models. The sum of all three 
positive Liapunov exponents is defined as "Kolmogorov entropy" 
(hk = In order to distinguish regular from stochastic or- 

bits for all models, unlike the case of Poon & Merritt (2004) who 
use the histogram of criTo, we consider the histogram of hkTr> 
for the orbits at each energy shell from both start-spaces. Since 
the histogram of h^To is found to have similar behaviours for all 
models, we present it in Fig. 1 for the orbits at shell 13 from the 
stationary start-space of the most stochastic model (i.e.. Models PT 
with 7 = 1.5). As can be seen from Fig. 1, the separation of orbits 
into two groups representing two peaks becomes apparent as the 
integration time increases. Out of these two peaks, one is situated 
at narrow regions near zero showing the regular orbits, whereas 
another one is located at non-zero value with larger spread that 
decreases with the integration time, representing the stochastic or- 
bits. After a time interval of ~ IO^Td, h^To of the stochastic orbits 
at each energy shell for all models are found to approach towards 
common values and their mean values do not change very much. 
This suggests the followings: (1) The stochastic orbits in all models 
diffuse slowly through their allowed phase-space, and can sustain 
their regular shapes for ~ IQ^To or longer. (2) The mixing-time 
that is associated with diffusion through the Arnold web would be 
~ KfTo (cf. Merritt & Valluri 1996), since it is clear from Fig. 1 
that there is still a lot of mixing going on at ~ IQ^To- Our com- 
putation of the Liapunov exponents for periods of IQr'To, which is 
two orders of magnitude longer than others in the hterature, allows 
us to achieve these results. 



At the integration time of IQ^Tu, the critical value, hicTo ~ 
10"'"^, is found to separate the orbits into two peaks in the his- 
tograms of hkTr, at each energy shell from both start-spaces for all 
models. Thus, the orbits are considered to be chaotic if HiJ'd > 
hkcTo- As can be seen from Fig. 2, a large fraction of the orbits 
at all energy shells from the stationary start-space are found to be 
stochastic for all models. For each of the two selected values of y, a 
larger fraction of the stochastic orbits are found in Models PT than 
in Models OT. Furthermore, Models PT (Models OT) with y = 1.5 
have a relatively larger fraction of the stochastic orbits than those 
of Models PT (Models OT) with y = 1.0. On the other hand, a 
significant fraction of the orbits from the x - z start-space at each 
energy shell are also found to be stochastic for all models. 



6. Construction of self-consistent models 

Self-consistent equilibrium models can be constructed by the usual 
method of Schwarzschild (1979). Here a linear superposition of or- 
bits is sought, each populated with an appropriate number of stars 
that could reproduce the mass of each cell. This method is formu- 
lated as 

M 

J^C(i)B(i,j) = D(j) (;=l...,Af) with 
(=1 

C(0>0, ii=l,...,M), (8) 

where B(i, j) is the time spent by the orbit in the cell, D(j) 
is the mass of the cell, and C(0 is a weight associated with the 
orbit, also called the non-negative occupation number of that 
orbit. To employ Schwarzschild's method, all models are divided 
into 960 cells by following the scheme set by Merritt & Fridman 
(1996). Furthermore, as in Schwarzschild (1993), the masses of 
960 cells and the time spent by each orbit in 960 cells are nor- 
malized to unity (i.e., li^^'J D(j) = 1 and ^^^l B(i,j) = 1). Using 
Lucy (1974)'s iterations method as formulated in Schwarzschild 
(1993), we have then minimized the mean square deviation in the 
cell masses, i.e., 

j=l V ,=1 / 

where = 960 and M - total number of supplied orbits. In this 
paper, we have used the departure from self-consistency to present 

our results, which is defined as 5 = ^ y[x^javerage mass per cellj 

(cf. Merritt & Fridman 1996). In order to start the Lucy iterations, 
we first choose the trial values of C(0=constant. Later, in each it- 
eration, we compute C{i)new to derive 6. We continue the iterations 
until satisfactory convergence in 6 is achieved, which is found to 
be around the ~ 30000 iterations for all models. We repeat this 
procedure for an ever-larger, randomly selected sample of orbits 
and then 5 as a function of the number of orbits is plotted. If a 
self-consistent solution exists, 6 would decrease rapidly with the 
number of orbits. 

7. Results and discussion 

Our result in Sect. 5 supports the findings of previous workers 
(Merritt & Fridman 1996; Wachlin & Ferraz-Mello 1998; Siopis 
& Kandrup 2000; Kandrup & Siopis 2003; Poon & Merritt 2004) 
that a large fraction of the orbits from the stationary start-space 
for models with y > 0.0 are stochastic, while there are signif- 
icant fraction of the stochastic orbits in the x - z start-space as 
well. Furthermore, a number of authors (Merritt & Fridman 1996; 
Holley-Bockelmann et al. 2001; Holley-Bockelmann et al. 2002; 
Poon & Merritt 2004) have claimed that cuspy triaxial equilibria 
can be constructed with sizable fractions of the stochastic orbits. 
This encouraged us to employ the stochastic orbits along with reg- 
ular ones while constructing the self-consistent solutions for all 
models Usted in Table 1. In order to carry out this, we follow the 
method given in Sect. 6, and the results are presented in Fig. 3. 
Since the mixing-time, ~ lO'^r^, is shorter than the integration 
time, IQ^To, we attempt to construct self-consistent solutions in 
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which the stochastic and regular orbits are treated in the same way, 
allowing each stochastic orbit to have an arbitrary occupation num- 
ber. Thus, a total of 4500 orbits are employed and then <5 as a func- 
tion of the number of orbits are plotted with "circles" in Fig. 3. 
Here 6 decreases very fast with the number of orbits and converges 
well for aU models. Therefore, we conclude that the self-consistent 
solutions are found to exist for all models considered in this paper. 

Furthermore, we have run an experiment to construct the self- 
consistent solutions utihzing only the regular orbits, although there 
is no obvious physical reason why nature would host only the reg- 
ular orbits. The results of this study are represented by "crosses" 
in Fig. 3. Except for Models OT with y = 1.0, the remaining 
three models do not show better convergence in S with increas- 
ing number of provided orbits. So, we conclude that the regular 
orbits can provide a sufficient variety of shapes to construct the 
self-consistent solution for Models OT with y = 1.0, which has a 
large fraction of regular orbits than other selected models (see Fig. 
2). 

Although we have already shown the existence of the self- 
consistent solutions by "circles" in Fig 3, it would still be interest- 
ing to construct the "fully-mixed" solutions for all models due to 
the reasons listed in Merritt & Fridman (1996). So, we, finally, at- 
tempt to construct the "fully-mixed" solutions by following Merritt 
& Fridman (1996). The results are shown in Fig. 3 by "square" 
and "star" when the stochastic orbits are "fully-mixed" at 10 and 
15 lowest energy shells, respectively. For each of these two cases, 
the convergence in 5 is found to be good, which suggests that the 
"fully-mixed" solutions can be constructed for all models. Thus, 
we conclude that the self-consistent solutions exist for all models 
considered in this paper. 
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Appendix A: The force derivatives 

Here we present the derivatives of the gravitational forces corresponding 
to the ZC96 models, which are calculated by the gradient of the three 
components of forces, i.e., F, with i = 1, 2, 3, in the Cartesian coordinates. 

The derivatives. Fa = can be written as 

/ X \^ / X \^ I 2x^ \ 

Fii = - uoD [-fj +udTt- l^^j Ti Ts - -- + T, T^\ 

+ ^ + — 7-3 + 3 - T2Tf, + ZTA—^ + T2T^ 
r r \r I \ r j 



I w, X 

+ 6 I ^ + — r4 



, for! = 1,2,3. 



(A.1) 



Furthermore, the cross derivatives, Fij = -g^ with / 9t j, are defined as 

_ XjXj [udd "d T'l „ „ - T'l 



Fij = 



r I r r- r 
3 r, /. T, 



+ 2 73 



3 72 / l2\ 

+ —^T(, + 3[2--|j n - 6T4 



for i&y= 1,2,3. (A.2) 



In Eqs. (A.1)-(A.2), the terms Ti, T2, T,, T4, Tj, Tg, and T, have fol- 
lowing forms: Ti = {ix] - x\- xfj, T2 = [x] - xfj, T3 = - -jf ), 

7-4 = - = (f? - ^ + 2^) = - + ^), and 
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Table 1. Model parameters and Liapunov time (t/,) 
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Model 


P 


1 


r 


Tl (in units of To) 


Comments 




0.65 


0.60 


1.0 


-21 ± 0.46 




Models PT 










Prolate Triaxial 




0.65 


0.60 


1.5 


-13 ±0.33 






0.95 


0.60 


1.0 


-25 ± 0.50 




Models OT 










Oblate Triaxial 




0.95 


0.60 


1.5 


-17 ± 0.40 





Table 2. Shell parameters for Models PT and OT 



Shell 


Radius " 

y=1.0 7=1.5 


Energy for Models PT 


Td for Models PT 


Energy for Models OT 


Td for Models OT 


7 = 1.0 


7 = 1.5 


7 = 1.0 


7= 1.5 


7 = 1.0 


7= 1.5 


7 = 1.0 


7 = 1.5 


1 


0.2791 


0.1512 


-0.8103 


-1.3203 


2.916 


1.059 


-0.7963 


-1.2974 


3.204 


1.156 


2 


0.4464 


0.2635 


-0.7265 


-1.1378 


3.980 


1.668 


-0.7089 


-1.1116 


4.356 


1.812 


3 


0.6076 


0.3760 


-0.6600 


-1.0073 


4.956 


2.260 


-0.6407 


-0.9800 


5.412 


2.452 


4 


0.7744 


0.4949 


-0.6023 


-0.9017 


5.964 


2.876 


-0.5824 


-0.8744 


6.492 


3.116 


5 


0.9530 


0.6238 


-0.5503 


-0.8112 


7.044 


3.548 


-0.5304 


-0.7845 


7.644 


3.836 


6 


1.1483 


0.7662 


-0.5024 


-0.7310 


8.244 


4.299 


-0.4830 


-0.7054 


8.916 


4.636 


7 


1.3660 


0.9259 


-0.4575 


-0.6584 


9.604 


5.156 


-0.4390 


-0.6343 


10.357 


5.548 


8 


1.6124 


1.1075 


-0.4151 


-0.5920 


11.173 


6.156 


-0.3978 


-0.5694 
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Here uod, vdo, and Wod are given by 



/ r 1 1 

\r+ 1/ \r+ I r 



+ (1-7) 



1 



vdd 



Wdd 



r, r 



(r + r^f-y 



r+ 1 

(l-7)(2-7)-2(2-7)(4-7) 



1 

r + \ r 
r 



+ (4 - 7)(5 - 7) 



r 



(l-7)(2-7)-2(2-7)(4-7) 



(r + r2) 



+ (4 - 7)(5 - 7) 



(r + r^r- 



(A.3) 



which represent the derivatives of wo, vo, and wo with respect to r, respec- 
tively. 
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Fig. 2. Fraction of the stochastic orbits per shell in the orbit li- 
braries. The left panels are for the stationary start-space, whereas 
the right ones are for the x - z start-space. The cusp parameter 7 is 
given in top left comer of each panel. Circles: Models PT; crosses: 
Models OT. 
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Fig. 3. Departure from self-consistency (6) as a function of the 
number of orbits supplied to the minimization routine. Left pan- 
els: Models PT with 7 - (1.0, 1.5); right panels: Models OT with 
7 = (1.0,1.5). Circles: all orbits; crosses: regular orbits only; 
squares: the stochastic orbits are fully mixed at 10 lowest energy 
shells; stars: the stochastic orbits are fully mixed at 15 lowest en- 
ergy shells. 



